Regional heatmaps of model performance under extremes

Overview

This notebook will guide you through the code to plot regional heatmaps of model performance (KGE and hit rate) under climatic extremes in R. The following figures will be computed:

• Figure 4: KGE regional heatmap for aggregated crops

• Supplementary figures 34-36: regional heatmaps of all KGE components for aggregated crops

• Supplementary figure 37: regional heatmap of the model hit rate for aggregated crops

• Supplementary figure 38: regions used in the heatmap indicated on a world map

The fully processed and filtered data files, available under GGCMI-validation/data/processed/figure_ready_data are required to produce the figures.

We load the necessary libraries

Code

library(tidyverse)   # Includes dplyr, ggplot2, readr, etc.
library(sjPlot)  # to set and configure plotting theme
library(rnaturalearth) # to create base world map with countries
library(spdep) # spatial functions
library(viridis) # color-blind proof color scheme

We configure the plotting theme based on the data visualization blog from Cédric Scherer: https://www.cedricscherer.com/

Code

theme_Scherer <- function (base_size = 12, base_family = "Helvetica") {
  half_line <- base_size/2
  theme(
    line = element_line(color = "black", linewidth = .5,
                        linetype = 1, lineend = "butt"),
    rect = element_rect(fill = "white", color = "black",
                        linewidth = .5, linetype = 1),
    text = element_text(family = base_family, face = "plain",
                        color = "black", size = base_size,
                        lineheight = .9, hjust = .5, vjust = .5,
                        angle = 0, margin = margin(), debug = FALSE),
    axis.line = element_blank(),
    axis.line.x = NULL,
    axis.line.y = NULL,
    axis.text = element_text(size = base_size * 1.1, color = "gray30"),
    axis.text.x = element_text(margin = margin(t = .8 * half_line/2),
                               vjust = 1),
    axis.text.x.top = element_text(margin = margin(b = .8 * half_line/2),
                                   vjust = 0),
    axis.text.y = element_text(margin = margin(r = .8 * half_line/2),
                               hjust = 1),
    axis.text.y.right = element_text(margin = margin(l = .8 * half_line/2),
                                     hjust = 0),
    axis.ticks = element_line(color = "gray30", linewidth = .7),
    axis.ticks.length = unit(half_line / 1.5, "pt"),
    axis.ticks.length.x = NULL,
    axis.ticks.length.x.top = NULL,
    axis.ticks.length.x.bottom = NULL,
    axis.ticks.length.y = NULL,
    axis.ticks.length.y.left = NULL,
    axis.ticks.length.y.right = NULL,
    axis.title.x = element_text(margin = margin(t = half_line),
                                vjust = 1, size = base_size * 1.3,
                                face = "bold"),
    axis.title.x.top = element_text(margin = margin(b = half_line),
                                    vjust = 0),
    axis.title.y = element_text(angle = 90, vjust = 1,
                                margin = margin(r = half_line),
                                size = base_size * 1.3, face = "bold"),
    axis.title.y.right = element_text(angle = -90, vjust = 0,
                                      margin = margin(l = half_line)),
    legend.background = element_rect(color = NA),
    legend.spacing = unit(.4, "cm"),
    legend.spacing.x = NULL,
    legend.spacing.y = NULL,
    legend.margin = margin(.2, .2, .2, .2, "cm"),
    legend.key = element_rect(fill = "gray95", color = "white"),
    legend.key.size = unit(1.2, "lines"),
    legend.key.height = NULL,
    legend.key.width = NULL,
    legend.text = element_text(size = rel(.8)),
    legend.text.align = NULL,
    legend.title = element_text(hjust = 0),
    legend.title.align = NULL,
    legend.position = "right",
    legend.direction = NULL,
    legend.justification = "center",
    legend.box = NULL,
    legend.box.margin = margin(0, 0, 0, 0, "cm"),
    legend.box.background = element_blank(),
    legend.box.spacing = unit(.4, "cm"),
    panel.background = element_rect(fill = "white", color = NA),
    panel.border = element_rect(color = "gray30",
                                fill = NA, linewidth = .7),
    panel.grid.major = element_line(color = "gray90", linewidth = 1),
    panel.grid.minor = element_line(color = "gray90", linewidth = .5,
                                    linetype = "dashed"),
    panel.spacing = unit(base_size, "pt"),
    panel.spacing.x = NULL,
    panel.spacing.y = NULL,
    panel.ontop = FALSE,
    strip.background = element_rect(fill = "white", color = "gray30"),
    strip.text = element_text(color = "black", size = base_size),
    strip.text.x = element_text(margin = margin(t = half_line,
                                                b = half_line)),
    strip.text.y = element_text(angle = -90,
                                margin = margin(l = half_line,
                                                r = half_line)),
    strip.text.y.left = element_text(angle = 90),
    strip.placement = "inside",
    strip.placement.x = NULL,
    strip.placement.y = NULL,
    strip.switch.pad.grid = unit(0.1, "cm"),
    strip.switch.pad.wrap = unit(0.1, "cm"),
    plot.background = element_rect(color = NA),
    plot.title = element_text(size = base_size * 1.8, hjust = .5,
                              vjust = 1, face = "bold",
                              margin = margin(b = half_line * 1.2)),
    plot.title.position = "panel",
    plot.subtitle = element_text(size = base_size * 1.3,
                                 hjust = .5, vjust = 1,
                                 margin = margin(b = half_line * .9)),
    plot.caption = element_text(size = rel(0.9), hjust = 1, vjust = 1,
                                margin = margin(t = half_line * .9)),
    plot.caption.position = "panel",
    plot.tag = element_text(size = rel(1.2), hjust = .5, vjust = .5),
    plot.tag.position = "topleft",
    plot.margin = margin(rep(base_size, 4)),
    complete = TRUE
  )
}
set_theme(theme_Scherer())

Then we can load our data:

Code

base_path <- "C:/Users/kobed/Documents/TOAD-PIKlegacy" # Where is the code repo stored?

crop_data <- readRDS(file.path(base_path, "extremes_aggr.rds"))

We have to join rnaturalearth regions to create these regional heatmaps

Code

# Load Natural Earth countries dataset
world <- ne_countries(scale = "medium", returnclass = "sf") %>% 
  mutate(ctr = iso_a3_eh) %>% # eh needed to get France and Norway as well 
  filter(name_long != " 
Ashmore and Cartier Islands", name_long != "Indian Ocean Territories") %>%  # Get rid of duplicate AUS entries
  dplyr::select(ctr,subregion) %>% 
  st_drop_geometry()  %>%  # Drop the geometry column because for AUS two different 
 distinct()  

crop_data <- crop_data %>% 
  mutate(ctr = ifelse(ctr == "JKX", "PAK", ctr)) %>%  # Replace "JKX" with "PAK" for Pakistan
  left_join(world, by = "ctr") %>% distinct()

Figure 4 - KGE heatmap

First we write a function to calculate the KGE. Note that we work with a customized version (KGE’) such that KGE = 1 - KGE’.

Code

# Function to calculate revised KGE components
calculate_kge<- function(benchmark, simulated) {
  r <- cor(benchmark, simulated, use = "complete.obs")  # Correlation
  Beta <- (mean(simulated, na.rm = TRUE) - mean(benchmark, na.rm = TRUE)) / sd(benchmark, na.rm = TRUE)  # Bias (Beta)
  Alpha <- sd(simulated, na.rm = TRUE) / sd(benchmark, na.rm = TRUE)  # Variability (Alpha)
  
  kge <- sqrt((r - 1)^2 + (Alpha - 1)^2 + (Beta)^2)  # Calculate KGE
  
  return(list(KGE = kge, r = r, Beta = Beta, Alpha = Alpha))  # Return all components
}

Then we calculate the KGE’ and its components for each climate extreme and model

Code

heatmap_data <- crop_data %>% 
  group_by(subregion, climate_extreme, model) %>% 
  summarise(
    KGE_components = calculate_kge(difftrend_obs, difftrend_sim)
  ) %>% 
  mutate(
    KGE = KGE_components$KGE,
    r = KGE_components$r,
    Beta = KGE_components$Beta,
    Alpha = KGE_components$Alpha
  ) %>% 
  dplyr::select(-KGE_components)  %>%  # Remove the list column, no longer needed
  distinct()

We define a model order to be visualised

Code

# Model order
order <- rev(c(
  "benchmark",
  "ensemble",
  "acea", 
  "simplace-lintul5",
  "dssat-pythia",
  "pepic",
  "cygma1p74",
   "ldndc",
  "pdssat",
  "epic-iiasa", 
  "lpj-guess", 
  "promet",
  "lpjml",
  "isam", 
  "crover"
))

Finally, we can visualise the heatmap. We use a logarithmic scale to better differentiate low and high values.

Code

heatmap_data$subregion <- factor(heatmap_data$subregion, 
                                 levels = rev(c("Western Europe", "Eastern Europe", "Southern Europe", "Northern Europe", "Australia and New Zealand", "Seven seas (open ocean)", "Northern America", "Central America", "Carribean", "South America", "Western Asia", "Central Asia", "Southern Asia", "Eastern Asia", "South-Eastern Asia", "Northern Africa", "Western Africa", "Middle Africa", "Eastern Africa", "Southern Africa"))  )

heatmap_data <- heatmap_data %>% # necessary because we call heatmap data later
na.omit()

heatmap_plot <- heatmap_data %>%
  filter(subregion != "Seven seas (open ocean)",
        subregion != "Carribean") %>% 
  na.omit() %>% 
   mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
  ggplot(aes(x = model, y = subregion, fill = -log(KGE))) + # We have to do - log(KGE) to get colorbar in right order
  geom_tile(color = "white", size = 1) + 
  
  # Add text annotations with the values
  geom_text(aes(label = round(KGE, 2)), color = "white", size = 3) +  # Round to 2 decimals
  
  # Set labels
  labs(x = "", y = "") +
  
  scale_fill_viridis(discrete = FALSE,  option = "viridis", direction = 1, name = "", 
                      breaks = c(-max(log(heatmap_data$KGE)), -min(log(heatmap_data$KGE))), # - to get order from low to high performance right
                     labels = c("Low Performance", "High Performance")
                    ) +
  
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1), # Rotate x-axis labels for readability
    axis.title.x = element_text(size = 12),
    axis.title.y = element_text(size = 12),
    legend.position = "bottom",  # Place the legend on the right for better presentation
    legend.key.width = unit(5, "cm"),  # Adjust the width of the legend key
    legend.key.height = unit(1, "cm"),  # Adjust the width of the legend key
    legend.text = element_text(size = 18),  # Set the size of the legend text
    legend.title = element_text(size = 10),  # Set the size of the legend title 
    legend.margin = margin(t = 10),  # Add some margin above the legend
    legend.box = "horizontal", # Ensure the legend spreads horizontally,
    strip.text = element_text(size = 18))  +          # Increase facet label size
  facet_wrap(~climate_extreme) 

heatmap_plot

Code

# Save the heatmap
ggsave(file.path(base_path, "fig5_HeatmapExtremes.png"), heatmap_plot, width = 15, height = 20, dpi = 300)

Supplementary figures 34-36: decomposition of the KGE heatmap

Alpha

On a log-scale

Code

heatmap_plot <- heatmap_data %>%
  filter(subregion != "Seven seas (open ocean)",
        subregion != "Carribean") %>% 
  na.omit() %>% 
   mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
  ggplot(aes(x = model, y = subregion, fill = log(Alpha))) + 
  geom_tile(color = "white", size = 1) + 
  
  # Add text annotations with the values
  geom_text(aes(label = round(KGE, 2)), color = "black", size = 3) +  # Round to 2 decimals
  
  # Set labels
  labs(x = "", y = "") +
  
  scale_fill_gradient2(
    high = "red",
    mid = "white",
    low = "blue",
    midpoint = 0,
    name = "",
    limits = c(min(log(heatmap_data$Alpha)), max(log(heatmap_data$Alpha))),
    breaks = c(min(log(heatmap_data$Alpha)), log(1), max(log(heatmap_data$Alpha))),
    labels = c("Lower variability", "Similar variability", "Higher variability")) +
  
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1), # Rotate x-axis labels for readability
    axis.title.x = element_text(size = 12),
    axis.title.y = element_text(size = 12),
    legend.position = "bottom",  # Place the legend on the right for better presentation
    legend.key.width = unit(5, "cm"),  # Adjust the width of the legend key
    legend.key.height = unit(1, "cm"),  # Adjust the width of the legend key
    legend.text = element_text(size = 18),  # Set the size of the legend text
    legend.title = element_text(size = 10),  # Set the size of the legend title 
    legend.margin = margin(t = 10),  # Add some margin above the legend
    legend.box = "horizontal", # Ensure the legend spreads horizontally,
    strip.text = element_text(size = 18))  +          # Increase facet label size
  facet_wrap(~climate_extreme) 

heatmap_plot

Code

# Save the heatmap
ggsave(file.path(base_path, "suppfig34_HeatmapExtremes.png"), heatmap_plot, width = 15, height = 20, dpi = 300)

Beta

Code

heatmap_plot <- heatmap_data %>%
  filter(subregion != "Seven seas (open ocean)",
        subregion != "Carribean") %>% 
  na.omit() %>% 
   mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
  ggplot(aes(x = model, y = subregion, fill = Beta)) + 
  geom_tile(color = "white", size = 1) + 
  
  # Add text annotations with the values
  geom_text(aes(label = round(Beta, 2)), color = "black", size = 3) +  # Round to 2 decimals
  
  # Set labels
  labs(x = "", y = "") +

   scale_fill_gradient2(
    high = "red",
    mid = "white",
    low = "blue",
    midpoint = 0,
    name = "",
    limits = c(-5 , 5),
    oob = scales::squish,  # squishes out-of-bounds values to limits
    breaks = c(-5, 0, 5), 
    labels = c("Overestimation bias", "No bias", "Underestimation bias")) +
    
  
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1), # Rotate x-axis labels for readability
    axis.title.x = element_text(size = 12),
    axis.title.y = element_text(size = 12),
    legend.position = "bottom",  # Place the legend on the right for better presentation
    legend.key.width = unit(5, "cm"),  # Adjust the width of the legend key
    legend.key.height = unit(1, "cm"),  # Adjust the width of the legend key
    legend.text = element_text(size = 18),  # Set the size of the legend text
    legend.title = element_text(size = 10),  # Set the size of the legend title 
    legend.margin = margin(t = 10),  # Add some margin above the legend
    legend.box = "horizontal", # Ensure the legend spreads horizontally,
    strip.text = element_text(size = 18))  +          # Increase facet label size
  facet_wrap(~climate_extreme) 

heatmap_plot

Code

# Save the heatmap
ggsave(file.path(base_path, "suppfig35_HeatmapExtremes.png"), heatmap_plot, width = 15, height = 20, dpi = 300)

r

Code

heatmap_plot <- heatmap_data %>%
  filter(subregion != "Seven seas (open ocean)",
        subregion != "Carribean") %>% 
  na.omit() %>% 
   mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
  ggplot(aes(x = model, y = subregion, fill = r)) + # We have to do - log(KGE) to get colorbar in right order
  geom_tile(color = "white", size = 1) + 
  
  # Add text annotations with the values
  geom_text(aes(label = round(r, 2)), color = "black", size = 3) +  # Round to 2 decimals
  
  # Set labels
  labs(x = "", y = "") +
  
   scale_fill_gradient2(
    high = "blue",
    mid = "white",
    low = "red",
    midpoint = 0,
    name = "",
    limits = c(-1 , 1),
    breaks = c(-1, 0, 1), 
    labels = c("Negatively correlated", "No correlation", "Positively correlated")) +
  
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1), # Rotate x-axis labels for readability
    axis.title.x = element_text(size = 12),
    axis.title.y = element_text(size = 12),
    legend.position = "bottom",  # Place the legend on the right for better presentation
    legend.key.width = unit(5, "cm"),  # Adjust the width of the legend key
    legend.key.height = unit(1, "cm"),  # Adjust the width of the legend key
    legend.text = element_text(size = 18),  # Set the size of the legend text
    legend.title = element_text(size = 10),  # Set the size of the legend title 
    legend.margin = margin(t = 10),  # Add some margin above the legend
    legend.box = "horizontal", # Ensure the legend spreads horizontally,
    strip.text = element_text(size = 18))  +          # Increase facet label size
  facet_wrap(~climate_extreme) 

heatmap_plot

Code

# Save the heatmap
ggsave(file.path(base_path, "suppfig36_HeatmapExtremes.png"), heatmap_plot, width = 15, height = 20, dpi = 300)

Supplementary figure 37: hit rate heatmap

First we define a function to calculate the hit rate

Code

calculate_hitrate<- function(benchmark, simulated) {
  
  # Check where the simulated values are less than 0
  hits <- simulated < 0
  
  # Calculate hit rate as the mean of hits
  hitrate <- mean(hits)  # This will give the proportion of TRUE values (simulated < 0)
  
  return(hitrate)  # Return all components
}

Then we can calculate per region and climate extreme

Code

heatmap_data <- crop_data %>% 
  group_by(subregion, climate_extreme, model) %>% 
  summarise(
    hitrate = calculate_hitrate(difftrend_obs, difftrend_sim)
  ) %>% 
  distinct()

`summarise()` has grouped output by 'subregion', 'climate_extreme'. You can
override using the `.groups` argument.

Now we can visualise the heatmap

Code

heatmap_data$subregion <- factor(heatmap_data$subregion, 
                                 levels = rev(c("Western Europe", "Eastern Europe", "Southern Europe", "Northern Europe", "Australia and New Zealand", "Seven seas (open ocean)", "Northern America", "Central America", "Carribean", "South America", "Western Asia", "Central Asia", "Southern Asia", "Eastern Asia", "South-Eastern Asia", "Northern Africa", "Western Africa", "Middle Africa", "Eastern Africa", "Southern Africa"))  )

heatmap_data <- heatmap_data %>% # necessary because we call heatmap data later
na.omit()

heatmap_plot <- heatmap_data %>%
  filter(subregion != "Seven seas (open ocean)",
        subregion != "Carribean") %>% 
  na.omit() %>% 
   mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
  ggplot(aes(x = model, y = subregion, fill = hitrate)) + # We have to do - log(KGE) to get colorbar in right order
  geom_tile(color = "white", size = 1) + 
  
  # Add text annotations with the values
  geom_text(aes(label = round(hitrate, 2)), color = "black", size = 3) +  # Round to 2 decimals
  
  # Set labels
  labs(x = "", y = "") +
  
  scale_fill_viridis(discrete = FALSE,  option = "viridis", direction = 1, name = "", 
                     limits = c(0,1),
                      breaks = c(0, 1),
                     labels = c("None captured", "All captured")
                    ) +
  
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1), # Rotate x-axis labels for readability
    axis.title.x = element_text(size = 12),
    axis.title.y = element_text(size = 12),
    legend.position = "bottom",  # Place the legend on the right for better presentation
    legend.key.width = unit(5, "cm"),  # Adjust the width of the legend key
    legend.key.height = unit(1, "cm"),  # Adjust the width of the legend key
    legend.text = element_text(size = 18),  # Set the size of the legend text
    legend.title = element_text(size = 10),  # Set the size of the legend title 
    legend.margin = margin(t = 10),  # Add some margin above the legend
    legend.box = "horizontal", # Ensure the legend spreads horizontally,
    strip.text = element_text(size = 18))  +          # Increase facet label size
  facet_wrap(~climate_extreme) 

heatmap_plot

Code

# Save the heatmap
ggsave(file.path(base_path, "suppfig37_HeatmapExtremes.png"), heatmap_plot, width = 15, height = 20, dpi = 300)

Supplementary figure 38: regions

Finally let’s visualise which grid cells belong to which regions from rnaturalearth

Code

world_geometries <- ne_countries(scale = "medium", returnclass = "sf") %>% 
  mutate(ctr = iso_a3_eh) %>% # eh needed to get France and Norway as well 
  filter(name_long != " 
Ashmore and Cartier Islands", name_long != "Indian Ocean Territories") %>%  # Get rid of duplicate AUS entries
  dplyr::select(ctr,subregion) %>% 
 distinct() 

regions_plot <- ggplot(world_geometries) +
  geom_sf(fill = "gray90", color = "black", size = 0.2) +  # Base world map
  geom_tile(data = crop_data, aes(x = lon, y = lat, fill = subregion), size = 1, alpha = 0.6) +  
  scale_fill_viridis_d(option = "C", name = "Subregion") +  # Improved color scale
labs(x = "", y = "") +
  theme(
    legend.position = "bottom",
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 10),
    plot.title = element_text(hjust = 0.5, size = 16)
  ) 

regions_plot

Code

# Save the plot 
ggsave(file.path(base_path, "suppfig38_Regions.png"), regions_plot, width = 28, height = 15, dpi = 300)